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Abstract 

An asymmetric stochastic process describing the avalanche dy- 
namics on a ring is proposed. A general kinetic equation which in- 
corporates the exclusion and avalanche processes is considered. The 
Bethe ansatz method is used to calculate the generating function 
for the total distance covered by all particles. It gives the average 
velocity of particles which exhibits a phase transition from an inter- 
mittent to continuous flow. We calculated also higher cumulants and 
the large deviation function for the particle flow. The latter has the 
universal form obtained earlier for the asymmetric exclusion process 
and conjectured to be common for all models of the Kardar-Parisi- 
Zhang universality class . 
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1 Introduction. 



Interacting particle systems with stochastic dynamics [1] and particularly 
the one-dimensional asymmetric exclusion process (ASEP) have been in- 
tensively studied [2]- [4], due to connections to growth processes [5], traf- 
fic flows [6], the noisy Burgers equation [7] and the Kardar-Parisi-Zhang 
(KPZ) equation [8]. Being one of the simplest examples of integrable non- 
equilibrium systems, the ASEP allows calculation of some dynamical prop- 
erties [9], [10], a large deviation function [11], and conditional probabilities 
[12]- 

In a standard formulation [13], particles move in such a way that there 
is always at most one particle per site on the one-dimensional lattice. Every 
particle hops to its right or left with biased probabilities provided the target 
site is empty. Otherwise, it does not move. Using a traffic terminology, 
this kind of interaction between particles can be called a "soft braking" . 

Another kind of interaction is an "aggressive braking" [14], [15], when 
a particle pushes the particle in front of it and then stops. The displaced 
particle shifts the next particle in front, if any, and so on. As a result, a 
chain of adjacent particles is shifted by one lattice space left or right at 
the same moment of time. Despite apparent non-locality of dynamics, the 
Bethe ansatz method is still applicable and the resulting Bcthc equations 
are solvable parallel to the ASEP case. Further generalizations of the ASEP 
have been proposed [16], [17]. In every case, however, an elementary motion 
of a particle produces a deterministic reconstruction (local or non-local) of 
the preceding lattice state. 

The beginning of intensive study of the ASEP nearly coincides with a 
burst of interest to the threshold dynamics and avalanche processes. Ap- 
peared originally in the sandpile model of self-organized criticality [18], 
the avalanche processes have been shown to be related to many different 
phenomena ranging from an interface dcpinning to earthquakes [19]. 

As an example of the threshold dynamics one can again consider ASEP- 
like stochastic model at the one dimensional lattice. In this case, however, 
we admit multiple occupation of a lattice site by particles. Like in the 
ASEP, each particle hops to its right or left. If the number of particles n 
at given site exceeds some critical value n c , the site is unstable and must 
relax immediately. The relaxation consists in spilling of m ^ n particles 
from the given site to neighboring sites by a fixed rule. If the neighboring 
sites become unstable, they relax as well. Thus, an avalanche of relaxations 
spreads over the lattice. The time interval between beginning and ending 
of every avalanche is negligible in comparison with characteristic hopping 
time of a single particle. 

Comparing to the ASEP with the aggressive braking, a fundamental dif- 
ference appears, when the spilling rule is stochastic [20]. In this case, the 
structure of avalanche becomes complicated. Unstable states may appear 
randomly even if an underlying structure of the lattice state is regular be- 
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fore an avalanche starts. Then, the distance at which avalanches propagate 
and the total mass of particles involved in an avalanche are random values 
described by probabilistic distributions. The configurations of particles in 
the lattice states before and after an avalanche may differ considerably, and 
the latter results from the first by a series of stochastic spillings. 

Another peculiarity of avalanche dynamics is a specific transition into 
a totally unstable state, when the density of particles exceeds some critical 
value and an avalanche never stops in the thermodynamic limit of infinitely 
large lattice [21]. This transition corresponds to change of the time scale 
characterizing the system, which can be defined for example as a ratio of 
system size to average velocity of particles. While for low density, the slow 
diffusion processes prevail, the fast avalanches dominate above the transi- 
tion point. The existing of two time scales was shown to be responsible 
for reaching of the self-organized critical state in systems with avalanche 
dynamics. 

The aim of this paper is to give a mathematical description of one kind 
of avalanche processes, the one-dimensional asymmetric avalanche process 
(ASAP) [22]. The ASAP is a partially asymmetric diffusion process with 
the totally asymmetric avalanche propagation. The similar directed one- 
dimensional stochastic avalanches have been considered in [23], where the 
asymptotic of avalanche distributions in the self-organized critical state 
have been calculated exactly for the open lattice in the thermodynamic 
limit. Instead, we study the ASAP on a ring with a fixed number of 
particles. In this case, the critical value of the density exists, depending 
on spilling probabilities, which corresponds to the transition from the in- 
termittent to continuous flow when the fast avalanche dynamics becomes 
dominating. In this paper, we concentrate on dynamical properties of the 
ASAP below this point. 

One of the reasons for the intensive interest to the ASEP is that it being 
exactly solvable gives a discrete version of the Kardar-Parisi-Zhang (KPZ) 
equation [10]. In the scaling limit, one can get analytically the universal 
quantities like critical exponents and scaling functions characterizing a vast 
class of nonequilibrium phenomena belonging to the KPZ universality class. 
On the other hand, the universal scaling properties of avalanche dynamics 
are much less investigated. There are very few successful attempts, [24]- 
[28], to find analytical arguments allowing one to relate the avalanche- like 
processes with one of well-defined universality classes such as the KPZ , 
Edwards- Wilkinson or directed percolation [29]. 

In the present work, we show that the generating function of the total 
distance Y t travelled by particles in the ASAP is given in the scaling limit 
by the expression 

l im \ / ^ lKl + K 2 G(K 3l ), (1.1) 

t^oo t 

where G{x) does not depend on parameters of the model and has the fol- 
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lowing parametric form 
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and Ki,K 2 , and K 3 are model dependent parameters. This universal form 
of function G(x) was claimed to be the feature of the KPZ universality 
class [30], [31]. It also determines the universal form of the large deviation 
function characterizing the deviations of the integrated particle current 
from its average value. These results give an evidence that the ASAP 
shares the KPZ universality class with the ASEP despite the significant 
difference in their dynamics. 

The article is organized as follows. In Section 2, we consider the master 
equation for a general stochastic model which leads to the ASEP and ASAP 
in particular cases. In Section 3, we derive the Bcthe ansatz equations 
for the generating function of the total displacement of all particles. A 
particular case of the ASAP is considered in Section 4 where the ASAP 
becomes identical to a drop-push version of the ASEP. The general case 
of the ASAP is considered in Section 5. Using the method developed in 
[32], [33], we analyse the integral equation corresponding to the density of 
roots of the Bethe equations and derive the generating function for the 
total displacement of particles. From it we obtain the expression for the 
cumulants of total distance travelled by particles, like mean velocity and 
variance, its large deviation function, and demonstrate that the ASAP 
belongs to the KPZ universality class. 3 

2 The master equation. 

In this section we are going to obtain the master equation describing ASAP, 
which is defined as follows. Consider the system of p particles on a ring of 
N sites as shown in Fig. 1. Particles jump left or right with probabilities 
Ldt or Rett, respectively, for infinitesimal time dt independently of each 
other. When a particle comes to already occupied site after the hopping 
either from left with the rate R or from right with the rate L , an avalanche 
starts. It develops step by step according to the following dynamical rules. 

If, at some step of the avalanche, n (n = 2,3, ...) particles are at site 
x, then 

with probability fj, n , n particles go to the site x + 1 ; 
with probability 1 — /i n , n — 1 particles go to the site x + 1 and one 
particle stays at the current site x. 

3 A brief summary of the results has been presented at " StatPhys-Taiwan 2002", 
during May 27 -June 1, see [34] 
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We imply that an avalanche takes infinitesimal time to end, i.e. from the 
point of view of Poissonian processes it plays a role of interaction resulting 
in the transition between configurations with single particle occupation. 
The totally asymmetric case discussed in [22] corresponds to particular 
choice of the rates, L = 0,R = 1. In the case of the ASEP a particle 
step to already occupied site is forbidden. However, it will be shown below 
to be closely connected with the ASAP. To make the presentation more 
systematic we start from the system of free particles then going to ASEP 
and ASAP dynamics. 

The state C of the system at time t is characterized by the probability 
Pt{C) satisfying the master equation 

d t P t {C) = Y,M{C,C')P t {C) (2.1) 

{C} 

The off-diagonal elements M(C, C) of the matrix M are rates of transi- 
tions from configurations C to C and therefore are always positive. The 
diagonal elements M(C, C) that give the total rate of the transition from 
the state C to all other configurations, enter the matrix M with a minus 
sign. Conservation of probability requires the identity 

M{C,C)=-^M{C',C) (2.2) 

{C} 

Let us return to particles at the lattice. Consider noninteracting par- 
ticles jumping left or right with probabilities Ldt or Rdt, respectively, for 
infinitesimal time dt . The probability Pt(x\, . . . , x p ) for particles to occupy 
sites x\, . . . ,x p obeys the master equation 

dtPt{x!, ...,x p ) = ~pP t (x 1} ...,x p )+ (2.3) 

p p 
Ly^Ptjxi, ...,Xi + l,...,x p ) + R.y^P t (xi, ...,Xj-l,...,x p ) 

if Xi + i — xi > l.We impose the condition L + R = 1 by an appropriate 
choice of time scale. 

In the ASEP, the form of the master equation should be modified if 
configuration C contains pairs of neighboring occupied sites. If there are 
only two neighboring sites, x, x+l, occupied by particles, the master equa- 
tion differs from Eq.(2.3) by the extra term LP t (. . . , x + 1, x + 1, . . .) + 
RPt{- ■ ■ , x, x, . . .) — P t (. . . , x, x + 1, . . .). If there are more than one pair, 
one must substract the other unwanted terms from the Eq.(2.3) for every 
pair to obtain equations taking into account the exclusion rules. Instead, 
one can reduce the consideration to the free master equation only, if one 
equates the appeared extra terms to zero putting the boundary conditions 
for the physical domain x\ < X2 < ... < x p : 

LP t {. . . , x+l, x+l, . . .)+RP t (. ..,x,x,...)-P t (...,x, x+l, . . .) = 0. (2.4) 
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The terms like Pt (. . . , x, x, . . .) corresponding to multiple occupation of 
sites do not contribute to the dynamics due to the exclusion rule. Therefore, 
they can be considered as auxiliary non-physical terms and redefined by 
the boundary conditions so that all extra terms in free equation vanish 
giving the correct equation for the system with interaction. 

Consider now more general condition of type Eq.(2.4) with the coeffi- 
cients a and (3 which do not coincide necessarily with the rates L and R in 
Eq.(2.3) and will be defined later, 

aP t (...,x,x,...)+pP t {...,x+l,x+l,...)-P t {...,x,x+l,...) =0 (2.5) 

To provide the probabilistic meaning of the Eq.(2.3) together with 
Eq.(2.5), some constraints should be imposed on a and (3 . In this case, 
we still use the exclusion rule that allows one to use the terms of type 
P t (. . . , x, x, . . .) as auxiliary ones which should be redefined in appropriate 
way. The condition Eq.(2.5) itself does not eliminate the contribution of 
extra terms yet. Nevertheless, we can try to use this condition to replace 
the unwanted terms by terms consisting of allowed configurations only. To 
this end, we can exploit the fact that two unphysical terms in Eq.(2.5) are 
of similar structure and consider this relation as a recursion: 

P t {...,x,x,...) = -P t (...,x,x + l,...)-Z-P t (...,x+l,x + l,...), (2.6) 
a a 

To proceed, it is convenient to consider the two-particle case separately. 
2.1 The case of two particles. 

If there are only two particles at the lattice, the recursion Eq.(2.6) can be 
immediately solved in terms of allowed configurations only: 

P t ( x , x ) = -J2(-^J P t (x-n-l,x-n) (2.7) 

where, due to periodic boundary conditions, all coordinates are integers 
modulo N . Substituting Eq.(2.7) into the Eq.(2.3) rewritten for the two 
particle case under the condition x 2 = X\ + 1, we get 

d t P t (x 1 x + 1) = LP t (x, x + 2) + RP t (x - l,x + 1) - P t (x, x + 1) (2.8) 




where [i = —a/ [3. To give a probabilistic meaning to the transition 
rates, the terms corresponding to the processes when the system leaves 
the configuration (x, x + 1) should be non-positive and those for coming 
into (x,x + 1) from other configurations should be non-negative. To keep 
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probabilities positive, we have to impose the condition that either \x is pos- 
itive or /j, = —R/L, when the term containing the infinite sum vanishes. In 
addition, conservation of probability, Eq.(2.2), requires 

a + (3= I. (2.9) 

Then, the condition n = —R/L implies a = R, [3 = L, i.e. the ordinary 
ASEP. In the case n > , we have 

a = -/i/(l-/i),0=l/(l-/i) ( 2 - 10 ) 
and can rewrite Eq.(2.8) in the form 

d t P t {x,x + 1) = LP t (x, x + 2) + RP t (x - l,x + 1) - P t (x,x + 1) (2.11) 

OO 

+ (R + L/j,)(-P t (x, x + 1) + (1 - n) ^ n ~ X Pt{x -n,x-n+l) 

n=l 

In terms of \x , the boundary condition, Eq.(2.6) reads 

P t (x,x) = (1 - n) P t (x - l,x) + nP t (x - l,x - 1) (2.12) 

The expression Eq.(2.11) shows that in addition to the Poissonian hopping 
given by the original kinetic equation, new terms appear in the equation 
which correspond to transitions to the configuration C = (x,x + 1) from 
the configurations {C} = {(x — n,x — n + 1), n = 1,2, ...}. The rates of 
the transitions are (R + L/i)(l — ii)/i n ~ 1 . In the case of two particles, these 
rates determine the avalanche dynamics defined above. 

To show the relation between ASEP and ASAP before going into details 
of the solution, let us consider the moment when an avalanche starts. This 
happens if either the left particle of two neighboring ones moves right with 
the rate R or the right particle moves to the left with the rate L and then 
two particles together make at least one step together with probability pi. 
Then an avalanche starts with the rate (R + Lp). This expression indeed 
enters those parts of kinetic equation, which correspond to the avalanche 
dynamics. Thus, the rate of beginning of an avalanche becomes zero when 
ji = —R/L, and only exclusion dynamics remains. Therefore, we may treat 
ASEP as an analytical continuation of ASAP with a parameter /i taking a 
special negative value. 

2.2 Many particle processes. 

One can expect that the n-particle interactions imposes n new constrains 
on the master equation Eq.(2.3). However, in this section we will show 
that under certain constraint on toppling probabilities /i„ no new bound- 
ary conditions appear and Eq.(2.12) is sufficient to take into account the 
interaction of arbitrary number of particles. 
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To generalize the boundary condition, Eq.(2.12), for the description of 
many particle dynamics defined above one should express an unstable con- 
figuration via configurations appeared at the previous steps of an avalanche. 
The form of these conditions depends on the fact whether the site x — 1 is 
occupied or not: 



P t (. ..,x- l, x, ,x , . . .) = (1 - Hn)Pt(- ■ ■ , x-l,...,x- 1 , . . .) + 

n— 1 n 

(1 - .., x-l,. x-l ,x,...) (2.13) 

n-1 

if the site x — 1 is occupied, and 

P t (... ,X, . ,X, ...) = HnPt{- ■ ■ , X~ 1,...,X- 1 , . . .) + 

n n 

/i„_iPt(...,a;-l,...,a:-l,a;,...) (2.14) 

V v ' 

n-1 

if x — 1 is empty. 

Like the two-particle boundary condition, Eq.(2.12), the many-particle 
conditions, Eqs.(2.13) and (2.14), should be applied recursively. Applying 
this recursion step by step to infinity, we generate an infinite series that 
consist of the transition probabilities between stable configurations only. 

On the other hand, one may treat the term P t (. . . , x, . . . , x, . . .) formally 

n 

applying the two-particle boundary condition Eq.(2.12), reducing sequen- 
tially the number of particles in unstable sites. As a result, we obtain 

P t (... ,x, ,x, . . .) = ii n P t {- ■ ■ , x- l,...,x- 1 , . . .) + 

n n 

(l-fi n )P t {...,x-l,...,x-l,x,...) (2.15) 

V v ' 

n-1 

Due to recursion, parameters \x n are expressed through the only parameter 
fi: 

M2 = /U, H3 = (l-fl)fl, /Un = (1 - V)Hn-l + W*n-2, n>3 (2.16) 

or 

M»=M \_r } ■ (2.17) 
1 + /i 

Generally, Eq.(2.15) and Eqs.(2.13),(2.14) do not coincide. However, the 
final series entering the kinetic equation which result from the sequential 
use of either latter or two former of them will be the same, provided that 
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the transition probabilities /i„ from Eqs.(2.13) and (2.14) governing the 
avalanche dynamics satisfy Eqs.(2.17) and (2.16). 

Indeed, every term of resulting series represents a finite avalanche. It 
contains the product of terms /i n and (1 — fi n ) coming from the successive 
use of the recurrent relations, Eqs.(2.13),(2.14) or Eq.(2.15). The first term 
of the r.h.s of Eq.(2.13) increases the number of particles at the unstable site 
by 1 in comparison with the l.h.s. The second term of the r.h.s of Eq.(2.14) 
decreases the number of particles by 1. Obviously, in every finite avalanche, 
the numbers of decreasing and increasing events are equal. Therefore, the 
coefficients (1 — /i„) from Eq.(2.13) and [i n -\ from Eq.(2.14) always enter 
the product corresponding to the avalanche in pairs and interchanging their 
places do not affect the structure of the final series. On the other hand, 
the interchanging /i n -i and (1 — ^„) between Eqs.(2.13) and (2.14) leads 
to Eq.(2.15) obtained from Eq.(2.12). 

Thus, we have shown that the two-particle boundary condition, Eq.(2.12), 
is sufficient for obtaining the kinetic equation for many particle avalanche 
process governed by the above dynamical rules. In the Bethe ansatz for- 
malism, a similar procedure is known as the two- particle reducibility [38] 
and provides integrability of a system. 

3 The Bethe Ansatz for generating function. 

Consider the generating function of the total path Y travelled by all par- 
ticles between time and t provided that the system is in configuration C 
at time t 

oc 

F t (C)= Pt{C,Y)e' Y , (3.1) 

Y=-oo 

where Pt(C, Y) is the joint probability that the system is in configuration 
C and the total distance is Y to the time t. The equation for the generating 
function and a boundary condition can be obtained from the master equa- 
tion Eq.(2.3) and the boundary condition for the probability Eq.(2.12) by 
multiplying by e 7 (e~ 7 ) every term corresponding to increasing (decreasing) 
distance Y by 1. 

p 

d t F t {xi, ...,x p )= Le^ ^2 Ft(xi, . . . , Xi + 1, . . . , x p ) + (3.2) 

i=l 

p 

+Re l ^ J F t (x 1 , ...,Xi-l,...,x p ) -pF t (x!, ...,x p ) 

i=l 

F t (...,x,x,...) = {l-fj,)e^F t (...,x-l,x,...)+ (3.3) 
He 2l F t {. . . ,x - 1, x - 1, . . .). 
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The sum of the function F t (C) over all configurations gives the generat- 
ing function of moments of total distance Y t , whose behavior is determined 
by the largest eigenvalue A(7) of Eq.(3.2) for large t . 



^F t (C) = (e^)~e A ^*. 

{C} 



(3.4) 



The derivatives of A(7) at 7 = give the cumulants of the distance Y t 

(Yt) 



lim 
lim 

t^oo t 



t^oo t 

(Y?) L 



lim ^ - ^ 

t— >oo t 



7=0 



lim 

t^oo t 



lim 

t— >oo 



lim 

t— >oo 



(if) - (Y t f a 2 A( 7 ) 



dj 2 



7=0 



<y t 3 ) + 2<y t ) 3 -3(y t 2 )<y t ) _ a 3 A( 7 ) 



Qj3 



7=0 



The quantity of our main interest is the large deviation function 



f(y) = um -InProb 

t^oo t 



(3.5) 
(3.6) 
(3.7) 

(3.8) 



which characterizes deviations of the distance Y t from the average value 
and can be expressed also trough the largest eigenvalue A (7) : 



m 

y 



(3.9) 
(3.10) 



Thus, we have to find the dependence of A (7) on the parameter 7. The 
master equation Eq.(3.2) and the boundary conditions Eq.(3.3) allow one 
to use the Bcthc ansatz in a usual form 



F t ( Xl ,...,x p ) = e At J2 A.{z ai ,...,z ap ) z ai xi 
CT (i p) 



(3.11) 



where the summation is over all permutations of (a\, . . . ,a p ) . The eigen- 
value corresponding to eigenvalue Eq.(3.11) is 



i=l i=l 

The parameters satisfy the Bethe equations 



z k 



i=i 



(1 - i^eJzj - fie^ZjZk 



(3.12) 



(3.13) 
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which follow from the substitution of Eq.(3.11) into Eq.(3.3) and the pe- 
riodic boundary conditions. The largest eigenvalue of the master equation 
corresponds to the stationary state of Markov process, so, one has to choose 
the solution of Eqs.(3.13) which provides the eigenvalue of the master equa- 
tion for probability of Eq.(2.3) to be equal to zero 

lim A( 7 ) = 0. (3.14) 

The Perron - Frobenius theorem ensures that this eigenvalue has no crossing 
in the whole range of 7. 

To solve the Bethe equations, it is convenient to transform the variables 
Zk- After the change of variables 

z k = ^—^e-i (3.15) 

the system (3.13) can be rewritten in the following way 



\ 1 + flXk/ 

Corresponding eigenvalue has the form 

1 - x k 



jJi x 3 + V x k 



^) = E(^ + ^)-, (3.7, 



fc=l 



4 The limit \l = 0. 

In the limit /j, — > 0, the ASAP becomes a particular case of the two- 
parameter family of exclusion processes discussed in [15] which, in turn, 
degenerates into the n = 1 drop-push model [14] in the totally asymmet- 
ric case L = 0. In this case, particles perform the partially asymmetric 
random walk with rates L and R. Going right, a particle jumps to the 
closest unoccupied site at its right overtaking all adjacent particles next to 
it, whereas the motion to left obeys the exclusion rule. 

The simple form of the Bethe equations in this case allows one to use 
the method proposed in [11] to obtain a full solution of the problem. If one 
introduces the parameter 

B = (-ir 1 e^ N f[x j (4.1) 

3=1 

the solution of Eqs.(3.16) will be given by the roots of the polynomial 
equation 

B (1 - x) N - x p = 0. (4.2) 
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To get the largest eigenvalue, one has to choose p roots approaching zero 
when 7^0. Following [11] one can use the Cauchy theorem to evaluate 
the sum over these roots integrating along the contour enclosing all the 
roots in a small vicinity of zero. 



i r -N- + E 



= -^/(.) (4.3) 

3— 1 x p 

Inserting R (1 — x) + L / (1 — x) — 1 instead of f(x), we get the expression 
for the eigenvalue in terms of series in powers of B 



AW = V B k (-l) kp C pfe - 1 ( - L< " 1 

hi Nk - 1 \k(l-p) + l/N k-l/N) 



(4.4) 



where p — p/N and C a — a\/(b\(a — 6)!) is the binomial coefficient. On 
the other hand, one can get the expression for 7 requiring z j = 1> 

which is correct for the groundstate solution . Taking the logarithm of this 
product and using Eq.(4.3) one gets 

P k=l 

Resolving two series of Eqs.(4.4), (4.5) and using Eqs.(3.5), (3.6) and (3.7), 
we obtain the expressions for cumulants of the total distance traveled by 
particles: 

R 



((l-p) + l/N) (2(l-p) + l/iV) 

2 



t^oo t 



... L-2&-1] 2 ( L(l-p) 



R 



((l-p) + l/N) (2(l-p) + l/iV) 

. cfST-i ( L(i- P ) 

[C P N --\] 3 V(l-V^)(3-l/iV) + 

R 

((l-p) + l/N) (3(l-p) + l/N) 

The scaling limit, N — > 00, of these expressions is of interest for us because 
it provides information about the large scale behavior independent of the 
details of microscopic dynamics. While the ASEP keeps the same universal 
behavior for any value of p, scaling properties of the ASAP may change 
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depending on how close to critical point the system is. In the case p = 0, 
the critical density corresponds to full occupation of the lattice, p c = 1. In 
the subcritical regime which corresponds to 

(1 - p) > 1/N (4.9) 

the situation is similar to the ASEP. In the subcritical region, the gen- 
erating function, A (7), takes the universal scaling form, Eq.(l.l) which 
has been already obtained for the ASEP and claimed to be universal for 
all models of KPZ universality class [11], [30]. Three model-dependent 
constants K\,K 2 ,K 3 in Eq.(l.l) are 

Kl = N "(jrh> <4 - 10a > 



K 3 = N 3/ V2tt (1 - p) p. (4.10c) 
The average velocity of particles, 

V x = - lim { ^^ W ?—-L{l-p), (4.11) 

p t^oo t (1 — p) 

and the other cumulants of the distance travelled by particles become diver- 
gent in the thermodynamic limit when p approaches 1. However, the physi- 
cal quantities characterizing the finite system should obviously be finite. As 
an example, one can consider simultaneous limit p — > 1,N — * 00. Substi- 
tuting, for instance, N^, (0 < i9 < 1) instead of (1-p) into Eqs. (4.10,4.11) 
one gets the expressions for cumulants of Y t which remain finite for finite 
N. The velocity of particles in this case, Voo ~ A* 5 , becomes explicitly 
dependent on N. Being finite for finite N, it is divergent when N tends to 
infinity. The upper limit for the exponent given by Eq.(4.9), is due to 
the term 1/N in Eq.(4.4). This term plays the role of "infrared cutoff" 
at the scale N, which ensures to remain of order N if (1 — p) becomes 
zero. When (1 — p) becomes of order of 1/N, the generating function, 
A (7), looses its universal structure, Eq.(l.l). Practically, this means that 
the presence of characteristic length breaks the scale invariance specific 
for the KPZ dynamics: 

N -» XN, K x - > \K U K 2 - > \~ 3/2 K 2 , K 3 - > X 3/2 K 3 , 

which is held in the subcritical region, Eq.(4.9). The character of particle 
motion near the critical line becomes strongly collective. Eventually, in the 
limit N = p, the process is equivalent to totally asymmetric diffusion of a 
single particle, with the distance, Y t , and time, t, rescaled as follows 

Y t Y t N, t -» tRN. 
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Although the solution for p = allows one to approach the vicinity of 
critical line, it seems to be very specific as it does not, in fact, involve the 
avalanche dynamics. However, as we will see below, it catches the basic 
universal scaling properties of the subcritical dynamics of the model for 
arbitrary p. 



5 The case of arbitrary ji < 1. 



In the case of arbitrary p, the Bethe equations Eq.(3.16) cannot be reduced 
to a polynomial equation. Nevertheless, they still can be solved in the limit 
N — > oo,p — > co, p = p/N = const. Let us consider the equation obtained 
by taking the logarithm of both parts of Eq.(3.16) 

1 P 

Po(xk) - j-f^2®( x ]/ x k) - 7 = 2iriZ(x k ) (5.1) 
i=i 

Q(y/x) = ]nx-]ny + ]n ]y y/ . X 

1 + fJ-x/y 

P0 = ln(7f^) (5.2) 
V 1 + px k ) 

We define <d (y/x) at the complex plane of variable y with branch cuts 
shown in Fig. (2). For small positive 7, the solution corresponding to the 
largest eigenvalue was shown in [30] to behave as 

Xfc^re 27 ™?, r ~ 7 1 ^ p , 7^0 (5-3) 

The radius r behaves nonanalytically when 7 approaches zero, so in the 
limit p — > 00 , radius r becomes finite no matter how 7 is close to zero. It 
can be easily verified that this is also correct for non-zero p at least in the 
limit 7 = 0. The analytical function Z(x) is fixed by the choice of logarithm 
branches. The distribution of roots given by Eq.(5.3) corresponds to the 
following choice: 

Then, assuming that the roots are at smooth contour at the complex plane, 
the derivative of Z(x) with minus sign has a meaning of density of roots 
along the contour 



Instead of Eq.(5.1), we are going to solve the equation 

p 

V 



1 p 

Po{x) - jz^Q{x /x) = 2iriZ{x) (5.6) 



i=i 
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together with Eqs. (5.4,5.5) under the assumption that the density of roots is 
an analytical function of x. To solve these equations one need to transform 
Eq.(5.6) to the integral form. This procedure is not straightforward and 
depends very much on properties of the function Z(x), which should be 
first assumed and then can be checked a posteriori. 

A simplest method is based on the replacement of the sum by the inte- 
gral along some contour in complex plane in the thermodynamic limit [35]. 
The analytic solution then can be found for some special cases, particularly 
for the case of the contour closed around zero [36] . It turns out to be the 
particular solution of our problem corresponding to single value of 7, 7 = 0. 
A development of this idea is the expansion method proposed in [37] for the 
investigation of the conical point of ferroelectric six-vertex model, which 
allows calculation of the leading term of the first cumulant of Y t [22] . To 
obtaine higher cumulants one needs to calculate the finite size corrections 
to the thermodynamic solution. This has been done for the ASEP [33] 
with the help of the method of perturbative expansion of the Bcthe equa- 
tions proposed by Kim [32] . For the ASAP we use the modification of his 
approach, which allows us to calculate the finite size corrections, avoiding 
some assumptions made in the original method. At least in the leading 
orders the results do not depend on these assumptions and we reproduce 
the results by Lee and Kim in the particular case \i = —R/L. To simplify 
the presentation we leave the details of the solution for the Appendices 
(A-C), going directly to the results. 

The solution of the Bethe equations results in the generating function 
obtained in the scaling limit 7A 3 / 2 = const 

A( 7 ) = 1 K 1 + K 2 G{ 1 K 3 ). (5.7) 

Here G(x) in a parametric form is defined by the relations 

G(x) = -Li 5/2 (-C), (5.8) 
x = Li 3/2 (-C), (5.9) 

and the function Life(x) is the polylogarithm defined by series 



n=l 

when |x| < 1. For arbitrary negative x , the integral definition can be used: 

00 



The equations, Eqs. (5. 7-5. 9), are nothing but Eqs. (1.1-1. 3) and Ki,K2, K3 
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'2tt 

(L-Ri-^y- 1 ) ( p V s \ s -l + 2p) 



are model dependent parameters 

K 2 = Ar- 3 / 2 i±^x (5.13) 

OO 

5 1 -ho' ((i - p)p) 3/2 

#3 = iV 3 /V2^(l - p)p. (5.14) 

The function G(x) has been already obtained for the ASEP [11] and 
claimed to be universal for all models of KPZ universality class [30]. Simple 
form of the eigenvalue (5.7) allows one to define a general expression for 
the cumulants of the integrated particle current in scaling limit 

limine = K 1 -K 2 K 3 =p(V av -V AS ep) (5.15) 

t— »oo t 

l im = K 2 KSG (n \0) = C<# + C<$ BP) n > 2, (5.16) 

t — »OC I 

where G( n \0) is n— th derivative of the function G(x) at x — 0. Here we 
divided the expressions for cumulants into two parts, which bring different 
physical content. Particularly Vav and Vasep, 



Lp + R 

(Pc - P) 



N (fie - p) 



(5.17) 



VASEP^L(l + p)(l-p)(l + j^j, (5.18) 

give the contributions to the average velocity coming from the avalanche 
part of the dynamics and its ASEP-like part respectively. The critical 
density p c , is defined as follows 

1 



1 + p 

which gives the phase diagram shown in Fig. 3, and the functions Fo (p, p c ) 
and F\ (p, p c ) arc nonsingular below and at the critical point, p = p Cl so that 
the form of Eq.(5.17) explicitly shows critical singularities of the average 
velocity. While Vasep gives the average velocity of the ASEP when p = 
—R/L, Vav vanishes at the same time. This is also the case for higher 
cumulants, which reproduce the results by Lee and Kim in this limit 

C<# ~ N 3J ^ L » + R (2np(l p))^ 1 pG™ (OjF^p.pe), (5.19) 

(Pc - P) 

C^sep ^(HA)^ (^(l-P))^ ) (0) . (5.20) 
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Thus, when the density of particles, p, approaches its critical value, p c , the 
divergency of average velocity of particles is characterized by the power law 

V~V AV ~ ( Pc - P y a (5.21) 

with the critical exponent a = 2 . In [22] this exponent has been shown 
to be nonuniversal with respect to the choice of different sets of dynam- 
ical rules p n . The other cumulants diverge as 4-th power of (p — p c ). 
The corrections to these laws are given in terms of the scaling variable, 
N _1 (p — pc) 2 ■ The condition this variable is small defines the limits of 
applicability of the perturbative scheme and bounds the subcritical region. 

\p-p c \ »i/Va. 

Closer to the critical line, we expect that the scaling will change as it did 
in the case p = 0. In addition, the case with arbitrary p allows one to 
consider the region of phase diagram above the critical line, where the 
average velocity of particles will grow with N. This case requiring the 
modification of finite size expansion scheme is not considered here. 

To obtain the cumulants, one needs knowing only the behavior of eigen- 
value in the vicinity of the point 7 = 0. At the same time, for the large 
deviation function, the whole range of 7 is relevant. The solution consid- 
ered above is valid for small negative 7. However, the definition of G(x), 
Eqs. (5.8,5.9), maintains analyticity in the region \C\ < 1. This implies 

Li 3/2 (-l) < 7*3 < Li3 /2 (1). (5.22) 

Beyond this domain one has to consider the solutions with different choice 
of the function Z(x). Instead, we can directly use the analytical continua- 
tion of G(x) proposed in [11]. To probe the whole range of negative 7, one 
can use the definition, Eqs. (5.8,5.9), where the functions Li 3 / 2 (C), Li 5 / 2 (C) 
are defined in the integral representation Eq.(5.11). For positive 7 outside 
of the domain Eq.(5.22), we use the following expression 

G(x) = ^[-ln(-C)] 3/2 -Li 5/2 (-C) (5.23) 

x = 4V^[-ln(-C)] 1/2 + Li 3/2 (-C), (5.24) 

where < C < — 1. Finally, using the definitions of G(x) in different 
domains of 7, we get the large deviation function in the scaling limit 

where the universal function H{x) is given by the following parametric 
expression 

H(x) = G(J3)-0G'(0) (5.26) 
x = G'(J3). (5.27) 
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In [11] this function has been shown to be skew, i.e. to have different 
asymptotic behavior for its argument taking on a large negative or pos- 
itive value. This is so for the ASEP, since the speeding up and slowing 
down are not equivalent due to the exclusion interaction. Specifically, it 
is much easier to slow down process by stopping or moving backward a 
single particle than to speed it up by moving forward all particles simul- 
taneously. Similar qualitative interpretation of this asymmetry for the 
ASAP is also possible. One can see from the explicit form of Eqs.(5.15- 
5.18), that the leading terms of contributions to average velocity coming 
from the ASEP and avalanche parts of dynamics have different signs, the 
positive term corresponding to avalanches and the negative being for the 
ASEP drift. This is why the universal function G (x) and subsequently 
H (x) are different from it's standard ASEP form in the minus sign before 
its arguments, Eqs.(1.2,1.3) (compare with Eqs. (20,21) from [11]). Thus, 
in our case the speeding up and slowing down are interchanged comparing 
to the ASEP. Indeed, to initiate an avalanche one particle should move 
faster or slower then the others to reach eventually an occupied site. To 
prevent an avalanche all particles should move simultaneously in the same 
direction. The latter has much less probability then the former. However, 
comparing to the ASEP the situation is even more peculiar. The solution 
of the ASEP is usually holds in the domain with definite direction of the 
drift, for example L < R. The opposite direction can be obtained by formal 
coordinate inversion x — > — x that is equivalent to L R. The solution 
of the ASAP holds for all values of L < 1 limited only by the condition 
\fi\ < 1. Thus, for small densities, p, and the rate of the left driving, L, 
close to 1 the situation may take place when Vasep > Vav, so that the 
resulting average velocity will be negative. This, however, does not affect 
the function H(x) as its skewness is related only to the avalanche direc- 
tion rather than to the Poissonian drift. At the same time full inversion 
transformation in the ASAP corresponds to simultaneous transformations 
L <-> R and fj, — ► l/fi, which changes the direction of avalanches as well as 
that of Poissonian drift. At the same time, if one considers the fluctuations 
of the particle flow in the reference frame attached to the average flow of 
particles the macroscopic fluctuations are quite similar to the ASEP up to 
change of the distance scale 



(Y t -(Y t )) 
(Y t -(Y t )) 



ASEP 



ASAP 



1 + 



/i + R/L 




(-tf s 2 (*-l + 2p) 
1-{-hY 2{{1- p)pf 



. (5.28) 



6 Summary and discussion. 
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To conclude, we have considered the asymmetric avalanche process on the 
ring. To introduce the avalanche dynamics to the master equation for 
the Poissonian process, we used the technique of the recurrent boundary 
conditions. We have solved the master equation by the Bethe ansatz and 
studied the solution corresponding to the stationary state. As a result, we 
have calculated the cumulants of the integral particle current, exactly in 
the case \i = and in the scaling limit for general \x. The large deviation 
function has been obtained, which has the structure, typical for models 
belonging to KPZ universality class. To calculate the finite size corrections, 
we used the modification of the perturbative scheme proposed by Kim in 
[32]. While in the leading orders the standard and modified approaches give 
the same results, it would be interesting to find out if this is so in arbitrary 
orders. Our investigation is valid for the densities below the critical point. 
However the model at the finite lattice can be considered for an arbitrary 
density of particles. One may expect that the scaling behavior of physical 
quantities should change above the critical point. The question of interest 
is how to modify the scheme to study the behavior of the model at the 
critical line and above. 
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A Solution of the Bethe equations. 

First, we suppose that the roots of the Bethe equation are placed along the 
closed contour T encircling zero, and the solution still preserves the invari- 
ance with respect to complex conjugation like the solution corresponding 
to 7 — > 0, Eq.(5.3). If we suppose also that the monovalucd analytical 
function Z^ 1 (?rj\r — -fa) can be defined everywhere at this contour, we ob- 
tain the mapping j — > xj, which allow us to use the Euler - Maclaurin 
summation formula for transformation of the sum over the roots of Bethe 
equations into an integral along the segment < j < p in the plane j with 
a correction term. The integral along the segment can be mapped into the 
integral along contour T in the plane x (see Appendix B) 

v v 

[f(Z- 1 (j))d ] + f. S .c^(ff(x)R(x)dx + f.s.c. (A.l) 

3=1 i 

where f.s.c. is the correction term. 

The functions of interest, e.g. the function 9 (y/x) can be represented 
on T as power series with additional logarithmic terms. In this case, the 
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only correction term appears which corresponds to a contribution originat- 
ing from the logarithm branch cut. As a result, the equation Eq.(5.6) 
simplifies to the following form 

X 

Po{x) - <j> 6 {y/x) R(y)dy - 7 = inp - 2ni J R{x)dx. (A.2) 

r x 

The reference point xo is the cross point of the contour with the posi- 
tive part of the real axis, which can be defined as x p e~ 2nl . The equation 
Eq.(A.2) coincides with that obtained in [22] in the limit of infinite lattice 
after replacing the sum by the integral and neglecting the finite size cor- 
rections. It turns out that the solution of Eq.(A.2) gives the exact solution 
of Bethe equations, provided that the inverse function Z _1 — is 
an analytical function in the segment 1 < j < p, what should be checked 
afterwards from the solution obtained. 

The only analytical solution of Eq.(A.2) corresponds to the case 7 = 0. 
Taking into account the normalization 



we get the solution 



R{x)dx = p. (A.3) 



2ni \x 1 — x 



According to Eq.(5.5) the function Z(x) corresponding to this solution, can 
be obtained by the integrating of the density R(x) : 

v I \ (P +1 ) 1 ( 1 x 1 1 - x 



27V 2m V x 1 - x 

Using the definition of Z(x), Eq.(5.4), we obtain the equation for the roots 
Xj of the Bethe equations. 

j = f,( p ln^-ln\^A (A.5) 
2m \ x 1-xqJ 

Considering j as a continuous parameter varying from to p, Eq. (A.5) 
can be treated as an implicit definition of the contour T. It is easy to check 
that Eq.(A.5) has a solution corresponding to a closed contour around zero 
when xq varies in the interval ^ xq ^ x p . The upper point x p , defined 
by the equation x p p — (1 — x p ) p p (l — p) 1_p = 0, is a monotonous function 
of p changing in the range ^ x p ^ 1/2 when p changes from to 1. 
Every value of xo less then x p corresponds to the particular position of the 
contour T passing through xo ■ As it was shown above, in the case of finite p 
and 7 = , all the roots are collapsed into one point x = 0. This situation 
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is realized when xq = 0. The nonzero values of xq correspond to two limits 
7 — > 0,p — > oo taken simultaneously. One can see from the Eq.(5.3) that 
if we put 7 ~ </> p , (1 < <f) < 0), the value of xo takes different limits when 
p — > oo depending on </>. At last, the case xq = x p is realized when the 
decay of 7 with growth of p is slower than exponential and particularly 
when the limit p — > 00 is taken for a fixed 7 and then 7 is put to zero. 
This case is of special interest for us because it gives the zero order of the 
solution we are looking for. When x > x p , the equation (A. 5) has no 
solutions corresponding to a closed contour encircling the zero point. 

Whereas the dependence of x p on p is defined by the irrational equation 
Eq.(A.5), for x = x p , the connection between p and the cross point of T 
with the negative part of real axes, x c is much simpler 

Importance of this point lies in vanishing of the root density Rq{x) at x c : 

R Q (x c ) = (A.7) 

This fact is crucial for the Kim's perturbative scheme [32] , [33] . Practically, 
it determines the range of applicability of Eq.(A.2). It has been noted 
above that Eq.(A.2) to be correct, the analytical function Z~ x (t^- — jj) 
should be defined to map the segment ^ j < p to contour T . It is possible 
only if the first derivative Z'(x) = — R(x) differs from zero everywhere in 
a given region [39]. According to Eq.(A.7), this is not the case at least at 
one point if 7 7^ 0. Thus, we cannot apply the above arguments to whole 
contour T. To overcome this difficulty, we can separate T into two parts 
and apply the summation formula to intervals where the requirements of 
analyticity are satisfied. It seems natural that all troubles with analyticity 
are concentrated around the point of crossing the contour with the negative 
part of the real axis. 

Consider the solution for odd p, and small positive 7. Then, the two 
roots closest to the negative part of the real axis, are conjugated to each 
other and can be denoted by 

x ( P -i)/2 = x c e~ le , X( p+ i)/2 = x c e le , (A.8) 

where x c and e are unknown real parameters. We exclude a small part of 
the contour T between Xr p _iy2 and £( p+1 )/2 an d assume the analyticity of 
Z- 1 (f^L _ i_) a t the segments 1 < j < (p - 1) /2 and (p + 1) /2 < j < p. 
This allows us to apply the Euler - Maclaurin formula for these segments 
separately. However, points j = (p + 1) /2 and j — (p — 1) /2 are located 
in the vicinity of a point on the real axis where Z~ x — jj) looses its 
analyticity. Therefore, the higher the order of derivatives of Z^ 1 — j^) 
at these points, the more singular their behavior. Practically, this means 
that all terms of the Euler-Maclaurin series have the same order in N. This 
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is why an alternative variant of the summation formula in the Abel-Plana 
form has been applied in [10]. The Abel-Plana formula requires, however, 
the analyticity of Z~ Y — j-f) in the strips 1 < Rej < (p — 1) /2 and 
(p + 1) /2 < Rej < p, which is assumed below. 

Using the Abel-Plana formula, we can transform the sum over the roots 
of Bethe equations into two integrals along the segments 1 < j < (p — 1) /2 
and (p + 1) /2 < j < p with a finite size correction term f.s.c. Then, these 
two integrals can be rewritten as the integral in the plane x along the closed 
contour T minus the integral along the small segment connecting the points 
x (p-i)/2 and iT(j, + i)/2. 



(P-1J/2 

- / / (z-H?")) # + / f (z-\j)) dj + f. s , 

1 (p+l)/2 
I 

j> f(x)R(x)dx- J f(y)R(y)dy + f.s.c. 

X(p-l)/2 

To get the expression for the finite size correction term, we have to know 
the behavior of the function Z^ 1 (£) in the vicinity of points £ = 0, 1/iV. 
Once we know Z(x), the inverse function can be constructed by inversion 
of its scries at every point. However, the inverse function Z -1 (£) become 
singular at the points £ = 0, 1/N in the limit N — > oo. This should be 
taken into account in constructing its series. 

Consider the Taylor series of Z(x) in points x — x c e Tie : 

oo ^ 

Z(z c e Tie + t) = y ^-t n 

We can introduce new shifted variables y = t — a T and consider the series 
expansion in y. 

Z(x c e^ -a T + y) = ^ + 5=F + £ (A.9) 



X (p+l)/2 



n! 

n=l 



where 



zZ = Z{x c e~ u ) = --; z + = Z(a: c e ie ) = 0; 

n=l fe=n v ' 

Two signs =F , marking all parameters here are to remind that we consider 
two expansions around points X( p _iy 2 = x c e^ te and X( p+1 )/ 2 = x c e le , and 
generally the parameters z~Z , <S~, 6^, cr_ are different from z^, <5+, b+, cr + . 
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The shift of parameters ct t is defined by the condition 

bf = (A.10) 

that determines the shift of the expansion to the point where the density 
R(x) is zero. Using the expansion (A. 9), we can construct the inverse series 



a = x c e+ + er, oi = W— , a 2 = -^79, 

03 " Ii7i(s)* (5!, »- 3!,A> '- 

Here, wc omitted the indices =p at variables a, 6, cr, <5 still implying two dif- 
ferent expansions. To ensure the proper choice of a branch of the square 
root in both expansions, one can check if Eq.(A.ll) for £ = is satis- 
fied. Finally after some algebra (see Appendix C) we rewrite the Bethc 
equations in the following form: 



R S = g^l)-l H*)"' (A.12) 

{ -L( T+ x° + -T-X°_)+R s ln?± + 

^ E — «~ " - *-"") + T^7 - T - x -) + 

2i^s-n K+ ' AN K + ; 



X 



if s 7^ and 



1 ~ / 1 \^ r(§ + i) 

AN^\2iN) 7rf +1 
[c- s (Li f +1 (-e— - ) - i"Li f + i(-e^-)) - 
c+ s (Li f+1 (-e— +) - (-if Li f+1 (- e -+))]} 



E„ „ sin es 1 , . s 

W— + p (A-13) 

s^0,s— — 00 
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otherwise. For 7 and p, we have respectively 



n=l v ' 

[c- (Li f+1 (-e— ") -rLi f+1 (-e--)) - 
c+(Lin +1 (- e — + )-(-i)"Lin +1 (- e - + ))] + 

(r_ lnx_ — r + lnx + ) + — (in 2 x_ — In 2 x + ) + 



■ — (t In T — Ti ltlli 1 J " ' 1 ' " •' " 

2N 

x+ n - xZ n x~ n \nx+ - xZ n \nx- 



EM 2 + 
^0 v n 



and 



p = 2niR . (A.15) 
Here we introduced notations t t , x T , s which are defined as follows 



.n=0 



^ T = -— (T^ + r T ); x T =x c e^ + ( TT (A . 16) 

OO / OO \ OO 

n=0 \n=0 ) n=0 



and R s are the expansion coefficients of Laurent series 



00 „ 

= E Jrr- ( A - 18 ) 



Using Eqs.(3.17),(B.7) together with Eqs.(C.l)-(C.4) we obtain the expres- 
sion of eigenvalue in terms of R s 

oc 

A( 7 ) = 2mN (-/*)"" A„i?„, (A.19) 

n=l 

where 

A„ = (i - i?^)"- 1 ) (1+m) 

While the equations Eqs.(A.12-A.15) look rather cumbersome, a significant 
simplification takes place if the following conditions are satisfied: 

7+ = r_ = r; (A.20a) 
x + = x- = x; (A. 20b) 

a~ = a+=a s . (A.20c) 

Then, instead of Eqs.(A.12,A.14), we get 
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R S = »(-«-D_J_ (-^'x 
2ni 2m i _ (- M ) |s| 

1 ^ / i Y T(n+ f) C2n+1,« T . , ^ 



n—0 x 



1 ^ / i \™ |) c 2 „+i . . . 



n— x / 

which together with Eq.(A.19) reproduces the results of paper [32]. The 
conditions of Eqs.(A.20a) and (A. 20b) are those accepted in [32] as assump- 
tions. Equation (A. 20a) is equivalent to equality S + = 5- (see Eq.(32) in 
[32]). Equation (A. 20b) is equivalent to the assumption that there is only 
one point x where Z'(x) = which is used as the expansion center in 
[32]. In our consideration, two complex conjugated points , x + and X-, 
are possible that merge into one point x in the limit p — > oo. The third 
equality Eq.(A.20c) is a direct consequence of first two. Generally, there 
are no obvious reasons for Eqs.(A.20a) and (A. 20b) to be satisfied .They 
can be checked a posteriori when the solution of Eqs.(A.12) and (A. 14) is 
obtained. We checked them in first orders of the perturbative solution and 
found that they are correct in the first three orders which are necessary to 
reproduce the results of [33]. 

To obtain the solution of the Eqs.(A.12) - (A. 15) which is consistent 
with the exact solution in the case /x = 0, one has to assume e to behave as 
e <~ N^ 1 / 2 when N — > oo. Therefore, we assume the following expansion 



(A.22) 



The other values in the Eqs.(A.12-A.15) can be represented as similar ex- 
pansions 

oo p(fc) oo oo 

k»* 

Equation (A. 4) is used as a zero order solution. Then Eqs.(A.12-A.15) 
should be solved order by order in powers of N~ x l 2 . The scaling depen- 
dence of 7 corresponds to 7./V 3 / 2 = const. The limit 7^0 corresponds 
to the limit ei — ► 0. The other parameters £2, £3, ■ ■ • depend on the way, 
how 7 approaches zero when N — > 00. However, the physical results do not 
depend on these parameters due to analyticity of eigenvalue. Solving the 
Eqs.(A.12-A.15) in first four orders, we get the expression for eigenvalue 
given in Eqs. (5.7-5.9) 
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B Evaluation of sums over roots of the Bethe 
equations. 

To evaluate the sum over roots of the Bethe equation, one can use the 
asymptotic formula approximating the sum by the integrals 

m 

m „ 

E /(•?') = / WW + corr U, n, m) (B.l) 
with the correction term given by the asymptotic Euler-Maclaurin series 

oo 

corr(f,n,m) = -(/(m) + /(„)) + ]T (/ ^ (m) - /^(n)) (B.2) 

or by the Abel-Plana integral form 

carr{f,n,m)= l -{f{m) + f{n))+ (B.3) 

OO 

1 f f(m + it) - /(m - it) - /(n + ft) + f(n - it) 

~ J ^T— x dt - (B-4) 

o 

The former requires the analyticity of the function /(j) at the segment of 
real axes j £ [n,m] and the latter at the strip of complex plane, Rex € 
[n, m]. Let us suppose, that the analytical structure of the function Z(x) 
allows one to define the analytical inverse function Z~ x — j^) that 
maps the segment j e [0,p] into the closed contour T encircling zero in 
the plane of the variable x. Then the derivatives with respect to j can be 
expressed in terms of Xj <XS follows 

i-J U) - -TT^TZT-J' (~ Z N + E i 1 ) ■ (B-5) 



8j J w NR{x 3l v . , 

We are interested in calculation of sums of the form 

where F(x) can be represented at the contour T as Laurent series with 
additional logarithmic term 

oo 

F(x) = Flnx + F ^ S - ( R? ) 

s= — oo 

The root x p lies at the real part of positive axes. We can introduce its 
"twin" at the other side of logarithm branch cut, Xq — G Xp , which cor- 
responds to j = 0. Then, application of Euler-Maclaurin formula Eq.(B.2) 
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gives 



JI E F ^ = | E - ) - / nv)R{y)dy + '^f (Bi 



Indeed, all derivatives in the Euler-Maclaurin series Eq.(B.2) taken at the 
points xq and x p , being equal,. are cancelled by each other. The only contri- 
bution to correction term comes from difference between imaginary parts of 
the logarithm at the banks of its branch cut. As it is shown in the section 
(5), this case is limited by 7 = 0. 

Let us consider the case when the roots 27j,_i)/2 and Xr p+ i\/i are lo- 
cated in the vicinity of the point x c satisfying Eq.(A.7), and, therefore, we 
can not guarantee existence of an analytical function Z^ 1 — jj) that 
maps the segment j € [0,p] to closed contour T. We, however, still assume, 
that the mapping like this exists at two its segments, which connect points 
xq, £(p_i)/2 and £( p +i)/i, x p . One can see from Eq.(B.5), that every deriva- 
tive with respect to j brings the competitive coefficients 1/N and 1/R(xj), 
i.e. when TV tends to infinity R (x( p ±i)/i) goes to zero. This is why we use 
Abel-Plana summation formula Eq.(B.3) instead of Eq.(B.2) to take into 
account all contributions of the same order in N. Applying it to each of 
two segments separately and using formula Eq.(A.8) for the roots £( p _i)/2 



and x 



(p+i)/ 



1 we obtain 



£ R S (I° + II) + ^F + 



(B.9) 



x° cos es x s 

— 1 

N N 



where 



2m ■ 



F Q + F(lnx +iir) s = ' 

p p x c sm(n-s)e , 

' ' : rn e(n-s) " 1 " 



-2ie < 



A ((In 



Ti 7^0, n=— 00 



sin ne 



(B.10) 



— coses) ,(B.ll) 
+ Flnx c s = 
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x s = 



z-H--) 



l + it 
N 



N 



N 

/(e™-l)dt, 



(B.12) 



oo 

i/{>^-'(-l^)-m Z -(-l^) 



N 



(B.13) 



-lnZ- 1 (-|)+lnZ- 1 (|)}/(e^-l)rfi 

and R s are the coefficients of the Laurent expansion of the density defined 
in Eq.(A.18). 



C Derivation of equations for R s . 

Rewriting the sum in the Eq.(5.6) with the help of Eq.(B.9) and collecting 
coefficient of the same powers of x we get 



s £ : 



2ni 



E 



n^s,n— — oo 



^i-(- M ) |s| 

— sin e(s — n) R n — cos es ■ 

s - n 2iN 



(C.l) 



2iN 



: 



eR = - J2 R 



s^O.s— — oo 

oo 

X, 



_ s sines 1 



7 = 2* ^ 

s^0,s — — oo 

p = 2iriRo 



sm es \ x 

e cos es i? 5 + — 

s \ s / N 



(C.2) 

(C.3) 
(C.4) 



where we use the function Q(y/x), treated as a function of the variable 
y, as the expansion valid at the contour T with coefficients defined like in 
Eq.(B.7) 

' s^O 



In a; s = 



9 = -1. 



27 



To go further one needs to obtain the explicit expressions for x s , x in terms 
of R s . Using the expansion of Z^ 1 + , Eq.(A.ll), we get 

*• = III (^Y KsK; - cl s Y+] (C.5) 

n=l ^ ' 

* = ^£(^v) 2 [U-W]> (C.6) 

n=l ^ ' 

where the coefficients and s are defined in Eq.(A.17) an the functions 
are given by 

_ i TjajjgT^gvwT^y^ 

z 7 e^* - 1 

o 

The method of evaluation of these integrals is described in detail in [33] . 
Finally, for we get 



r(| + i) 

7rf +1 



(Lin +1 (- e — t) _ (±») n Li» +1 (-e'"*))} . 



Using the equalities 



71=1 



Y^cl s {±i^^i) n { — \ =[Z-\4)]'-<$ t „ (C.9) 
f>T {±iy/^if (-^) 2 =\n[Z-\-4)\ -cj, (CIO) 

E OV^)" (-U * ¥^ = fn<$, + (en) 

n— 1 ^ / 2 

TV y x s R{x)dx, 



°0,s 



XX (±*V^)" (a^)' ?Jr = + (c - 12) 

n— 1 ^ '2 

TV y lna;i?(a;)rfa 



°0,s 
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Z-\z^) = x c e^ u ; cl s =x%; (C.13) 

and Eqs.(A.16) we come to the system of equations for R s , Eqs.(A.12)- 
(A.14). 
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Figure 1: The asymmetric avalanche process. 



Figure 2: The analytical structure of the contour T. Zigzag lines show the 
branch cuts of the function <d(y/x) in the complex plane of the variable 
y. The broken segment of the contour should be excluded from integration 
when 7 7^ 0. 



Figure 3: The phase diagram of the asymmetric avalanche process. 
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Fig. 2 
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